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Abstract Within pNRQCD we compute the masses 
of spin-averaged triply heavy baryons using the now- 
available NNLO pNRQCD potentials and three-body 
variational approach. We focus in particular on the role 
of the purely three-body interaction in perturbation 
theory. This we find to be reasonably small and of the 
order 25 MeV 

Our prediction for the fl ccc baryon mass is 4900(250) 
in keeping with other approaches. We propose to search 
for this hitherto unobserved state at B factories by ex- 
amining the end point of the recoil spectrum against 
triple charm. 

PACS 14.20.Mr ■ 14.20.Lq ■ 12.38.Bx 



1 Introduction 

Amongst the staples of hadron physics is Baryon spec- 
troscopy. Here, quark model computations of the light 
baryon spectrum [1,2,3] find only mild success beyond 
ground-states in various channels due to the plethora 
of open thresholds and couplings between channels. A 
much cleaner system can be provided by baryons com- 
posed of three heavy quarks (i.e. the combinations ccc, 
bbb, ccb and bbc). Since a considerable number of states 
are supposed to be below any strong decay thresholds, 
one can straightforwardly apply few-body reasoning, 
and quark model techniques can handle bound states 
better. 

Indeed, a review of the literature reveals many stud- 
ies that have computed or constrained heavy baryon 
masses, particularly the ground state fl ccc . We collect 
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some of the values obtained in section 7 below. Like its 
meson (quarkonium) counterpart [4,5], we expect this 
triply heavy baryon to attract much interest. 

For heavy quark systems, the development of poten- 
tial Non-Relativistic Quantum Chromodynamics (pN- 
RQCD) as an effective theory of QCD has allowed a 
more systematic treatment of the theoretical uncertain- 
ties involved in spectroscopic predictions by expand- 
ing in powers of 1/m [6,7,8]. For those in their ground 
state, pNRQCD can additionally be organized in stan- 
dard perturbation theory as a power expansion in a s [6, 
9]. While the theory itself has limitations due to the 
finiteness of the quark masses, the two-body static po- 
tential (quickly reviewed in section 2) has shown to be 
a good starting point for many meson investigations. 

With the three-body potential in NNLO perturba- 
tion theory now at hand [10] (we will give it in sec- 
tion 3), it is timely to perform an exploratory study of 
the ground-state triply heavy quark spectrum. This we 
present in section 8. The necessary QCD parameters a s , 
m c , nib are fixed by describing several common meson 
spectroscopy observables as explained in section 5. 

Finally, we comment on the feasibility of detecting 
the f2 ccc in section 10. Some numerical methods are 
relegated to the appendix. 



2 Static quarkonium potential in pNRQCD 

The static two-body potential for bound states of a 
quark and anti-quark is well known to NNLO (and be- 
yond) [11,12,13] as 
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The leading order potential is just the color Coulomb 
potential 



(o)_ 4 a s (r~ 2 ) 
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whereas the NLO and NNLO are, respectively, 
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By standard convention 7^ ~ 0.57721 ... is the Eulcr- 
Mascheroni constant; for three colors and in terms of 
the number of quark flavors Nf below the renormaliza- 
tion scale, the beta function that determines the run- 
ning of the coupling is expanded as 
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and the coefficients in the potential that remain in a 
conformal theory are 
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where £(3) ~ 1.202 ... is Riemann's zeta function. 

The NLO potential can be understood not only from 
pNRQCD (first, construct an effective theory around 
the heavy quark limit, then use Coulomb gauge in in- 
termediate steps to obtain the gauge-invariant poten- 
tial that is a matching coefficient in pNRQCD) but also 
from the D CTO -(|q|) time-like gluon propagator obtained 
in Coulomb gauge by Watson and Reinhardt [14], re- 
ducing it to its simpler Heavy Quark limit (there, sev- 
eral terms do not contribute to the Wilson loop poten- 
tial) 1 . 

In momentum space, the potential can be given to 
order NNLO [ ] at an arbitrary renormalization scale 
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with 

Cf^{x)=a 1 +^\og{x) , 

Cj^(x) = a 2 + Pi log 2 (x) + + 2A)Oi) Iog(x) . 

This obviously simplifies if the renormalization scale is 
chosen as q 2 itself 



V(q 2 ) = 
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2.1 1/m Potential 

In QED the 1 jm corrections to the static potential were 
long ago shown to vanish 2 . 

In perturbative QCD the non-Abelian vertex correc- 
tion produces a 1/m potential that cannot be gauged 
away. It can be nominally assigned to the 1/m 2 order 
through a field redefinition [7]. Since a recent study of 
the meson spectrum [16] finds reasonably large effects 
for the 15* states, especially in charmonium, we also 
consider it here. 

At leading order, the 1/m potential vanishes. For 
NLO and NNLO we employ the convention of [17]. Al- 
ternatively, one can use the NLO result (see Eq. (11) 
below) without [7] the factor (7/9), in order to match 
a lattice computation. 

In coordinate space the potential reads 



V m -i 



(11) 
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where m r is the reduced (pole) mass of the qq system, 
62 — —20.836 for Nf = 3 (appropriate for charmonium) 
and 62 — —18.943 for Nf = 4 (appropriate for bottomo- 
nium) and b 2 — —17.049 for higher scales where Nf = 5 
are given in [17]. The last term with a logarithm van- 
ishes if the scale is chosen as the BLM scale defined by 
Eq. (49) in the appendix. We have performed compu- 
tations with both this running scale and a fixed scale 
(m 2 or m 2 ). 

For example, in [15] it is shown that the Bethe-Salpeter 
ladder approximation generates Feynman gauge 1/m terms 
that vanish upon including the crossed-ladder box. No such 
terms are present in Coulomb gauge. 
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If the potential is constructed in momentum space, 
the 1 jm correction to the central static potential reads 



V m -i 
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Again, judicious choice of the scale fx = q disposes of 
the logarithm. 

A counterintuitive result is that matrix elements of 
the V m -i potential can actually be similar or, in ex- 
treme cases, even larger for bottom systems than for 
charm systems, since in a Coulombic system all ener- 
gies scale with the reduced mass m r . 

To see this, let us restrict ourselves to NLO and 
employ the convention of [7] in momentum space 



-2tt 



m r q 
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and compute (i/j\V m -i with a hydrogen like wave- 
function. Taking the Fourier transform of a IS" state 
with Bohr radius a 1 = m r a s yields 

1>{q) 
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and therefore 
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Extracting the dimensions and coupling, substituting 
the Bohr radius, and not minding about two constant 
positive numerical coefficients c, c', we find 



(Kn- 1 ) = —m r a A s {a s m r c)c' 



(16) 



If the coupling constant did not run, the expectation 
value (though suppressed in the perturbative counting) 
would be some factor of mb/m c ~ 3 larger for bottom 
than for charm systems. The NLO running of a s how- 
ever tames this growth and we will find very modest 
increases as a function of quark mass between charm 
and bottom. 

To give an example let's take pole quark masses of 
m c = 1.95 GeV, m b = 5.14 GeV and a s (m z ) = 0.114 
(various fits of these quantities will be given later on in 
section 6, these serve as illustration). 

At the soft scale, with reduced quark mass solv- 
ing iteratively as shown in appendix B.3, a s (^cts) = 
0.625 and a s ('-fas) = 0.395 Then a*(c)/a*(b) ~ 6.3. 
If the constant c in Eq. (16) is somewhat smaller than 
one, this number could be smaller and around 3. 



2.2 Running of the strong coupling constant 

The renormalization group equation that determines 
the running of the strong coupling constant to NNLO 
is 



dc 



d log q 2 
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By keeping the first term on the right hand side, or both 
terms, this equation can be solved to NLO or NNLO 
respectively. 

In general we employ the Runge-Kutta algorithm to 
numerically solve Eq. (17). To NLO the equation is also 
very simply analytically solvable and provides a handy 
check for the computer programme. Following [18] we 
introduce a scale A as is customary, so that 



NLO(r\1 
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with 



fa = 33 - 2N f 
4tt 12tt 

Inverting the equation yields 



A 2 = Q 2 exp - 



I 



ba* LO {Q 2 ) 



(19) 



In Table 1 we give for convenience, and as benchmarks, 
the values obtained by running back to low scales the 
renormalization group equation from the Z-boson pole, 
where the coupling constant is very accurately con- 
strained by many analyses [19]. At each of the scales 
mentioned in the table the number of active flavors in 
the beta function is decreased by one in a stepwise fash- 
ion and continuous matching is performed 3 . A typical 
run in agreement with world average and low-scale r 
data is plotted in figure 1. In addition we employ a 
second coupling that provides a best fit to key charmo- 
nium and bottomonium data, that while still broadly 
consistent with high energy data, is somewhat smaller. 

Moreover, we vary the number of flavors dynami- 
cally in the computer programme upon crossing each 
threshold. Since the quark masses themselves are vary- 
ing during each fit we cannot quote these thresholds 
here. However, they are in the vicinity of 1.6/1.7 GeV 
for the charm and 5 GeV for the bottom thresholds. 

We know that the value of a s at the Z-boson pole, 
evolved by backward NNLO running, is consistent with 

This introduces a non-analyticity that can, at least at 
NLO, be avoided by the use of the Brodsky-Lepage-Mackenzie 
method. See appendix A. This non-analyticity is enhanced if 
one employs the discontinuous matching conditions [20] based 
on effective theory, that we also intend to incorporate in fu- 
ture work. 
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Table 1 Benchmarks for a s (fi 2 ) at various scales. Input is 
the value of the coupling at the Z mass [19]. Shown are NLO 
and NNLO results with the number of active flavors decreased 
in a step-wise fashion from five to four and then three at each 
of the benchmark thresholds.. 



ft N f 


NLO NNLO 


91.188 5 
5 4 
1.6 3 
0.8 


0.1184(7) 0.1184(7) 
0.204(2) 0.2136(6) 
0.295(4) 0.336(2) 
0.417(8) 0.574(7) 




o.i- 



nuj 

0.1 1 10 100 1000 

q 2 (GeV 2 ) 

Fig. 1 Typical running coupling constant to order NNLO 
agreeing with recent r data and the world average [19]. 



the very precise determination based on radiative de- 
cays of the T meson [2 I], 0.184±0.015 at the bottomo- 
nium scale. Therefore, deviations from the Z-pole value 
in our fit results give us an idea of the size of the errors 
associated with the heavy-quark effective potentials. No 
further perturbation-theory approximation is implicit, 
given that the Schrocdinger equation is exactly diago- 
nalized. 

Finally, we remark that infinitely heavy quarkonium 
should be less sensitive to the infrared details of the 
interaction, but that some sizeable sensitivity remains 
because the quark mass is finite. Thus we employ the 
heavy quark-potential for all r in the programme. To 
avoid encountering spurious Landau pole singularities, 
we freeze the coupling constant at a low scale (400- 
600 MeV) and check the sensitivity to this procedure 
below. 

In conclusion, each of the runs reported will employ 
a slightly different fit value of a s , but typical shapes for 
the coupling constant can be seen in figures 1 and 2. 




0^ — — — — — 

0.1 1 10 100 1000 

q 2 (GeV 2 ) 

Fig. 2 Typical running coupling fit to two-body data. The 
actual best fit scale will depend on the order of perturbation 
theory, possible renormalon subtraction and treatment of the 
infrared, whether the static potential is or is not corrected by 
the 1/m force, etc., and varies from computation to compu- 
tation as will be indicated below. 

3 Heavy baryon potential in perturbation 
theory 

The static potential between three heavy quarks of equal 
mass in positions ri, r2, r3 is expanded in powers of the 
strong coupling constant as 

^(n.ra.ra) - V$ + + V<% LO + ■■■ . (20) 

In this work we will obtain the masses of a few heavy 
baryons from the Leading Order and Next to Leading 
Order potentials by employing a variational basis. Then 
we will study the effect of the intrinsic three-body force 
(star-shaped) that appears first at NNLO, and finally 
employ the rest of the known two-body NNLO terms 
to estimate further corrections to the mass values. 

The Leading Order potential is Z\-shaped i.e. given 
by the sum of the two-body Coulomb interactions 

LO 3 V|rr-r 2 | |r 2 -r 8 | + \r a -Ti\J " 

(21) 

For the remainder of this section we will shorten the 
notation by summing over an index i = 1,2,3 that runs 
over the three possible pairings of the quarks, so that 

i 

The one-loop corrections to this potential yield the 
NLO part. The coupling constant is renormalizcd and 
one needs to choose the renormalization scale at which 
the constant is initially given. Following the pNRQCD 
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Fig. 3 One of twelve Feynman diagrams contributing to the 
intrinsic three-body force in a triply heavy baryon at NNLO. 



custom, we first select the renormalization scale [x\ 
\j\vi | 2 . Then the potential to NLO [10] reads 



v LO + v NLO 



(23) 
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In momentum space, the NLO potential is easily 
reconstructed by comparison with Eq. (10). 

We now turn to the potential at NNLO. We are first 
of all interested in the intrinsic three-body piece, the 
star-shaped part of the potential that appears at this 
order. While three-body forces have been considered in 
the context of heavy hybrid mesons [22], applications 
to triply heavy baryons are sparse. 

This three-body force is conveniently organized in 
terms of an auxiliary potential 



V, 



NNLO-3 = 2(Kux(r 2 ,r 3 ) + V aux (ri,-r 3 ) 



+ V aux (-r 2 ,-r 1 )) , (24) 

that is computed via the Fourier transform 
dq 2 cfq 3 



Va UX (r 2 ,r 3 ) = i 
of a potential 

t-aux (02,03) = 
|q 2 + q3 



(2tt) 6 
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The two momenta g 2 and g 3 are flowing out of the two 
quark lines. 



The 2-body (zi-likc) contribution to the potential 
at NNLO is a simple generalization of Eq. (4) 



V, 
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That ground state heavy baryons feel more of a two- 
body like Z\-shaped rather than the T-shaped potential 
is supported by lattice data, where the A ansatz seems 
dominant up to distances of R ~ 0.7 fm [23] or even 
1 fm [24]. 

4 Infrared subtracted potential 

The static potential in terms of the pole mass is afflicted 
by an infrared singularity in perturbation theory [25, 
26,27] that can be subtracted by a momentum-space 
cutoff [28]. In passing from momentum space to coor- 
dinate space this amounts to a restricted Fourier trans- 
form 

d 3 q 

(27T) 3 ' 



V PS (r,Hf) 



r niol) 



(29) 



'|q|>M/ 

Correspondingly, the quark mass one deduces from the 
fits to data is in this particular "Potential Subtracted" 
scheme (shortened "PS" in what follows). The infrared 
singularity (and accordingly, bad behavior of perturba- 
tion theory) is now exposed as a counterterm necessary 
if one wants to retrieve the pole mass 

d 3 q 



:V{H) 



(30) 



|q|<« ( 2 ^) 3 

In perturbation theory, lacking another scale except the 
renormalization scale, the counterterm has to be pro- 
portional to [if up to logarithms, and is displayed ex- 
plicitly in Eq. (34) below. This is so that the static 
energy of the system does not depend on the subtrac- 
tion, 

V PS + 2m PS = V + 2m . (31) 

One can avoid the singularity altogether by directly re- 
lating the PS mass with the MS mass by means of 
Eqs. (32) and (34). 

The relation between pole quark masses m and MS 
masses m to NNLO is given in terms of the number of 
light quarks Nf as 

m = m(m) x (32) 



1 



4 a s (m) 



(13.44 



1.041JV/) ( 
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where a s is understood as the MS coupling constant 
(in terms of which we have anyway expressed the heavy 
quark static potential). The inverse relation is [29] 



m(m) = tox 

4 a s (to) 



(33) 



-14.33+ 1.0417V/) 



a s (to) x 



While using the Potential Subtraction scheme we 
need to convert the PS mass mps to the pole mass too. 
This is achieved by means of [28] 



to =mps(Hf) + [i/Cf x 



(34) 



4tt 
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5 Meson observables 

If we handle the heavy charmonium and heavy bot- 
tomonium systems we need to fix three parameters, 
the strong coupling constant at one scale a s (/x 2 ) (the 
evolution equation of the renormalization group then 
provides it at any other scale) and the quark masses 
at their pole (or in any other scheme). Alternatively, 
we employ a fixed running coupling constant from fig- 
ure 1 and fit only the quark masses to the ground state 
quarkonia (thus limiting ourselves to two parameters). 

We will employ five observables to fix these param- 
eters, two from charmonium (15 mass and radiative 
transition J/ip — > ij c j), two from bottomonium (IS 
and IP masses), and the B c mass. This will over con- 
strain the parameters for additional security and allow 
us to check the running of the coupling between the two 
scales, at the charm and at the bottom. 

Since the purely static potential docs not accommo- 
date a hypcrfine SrS2, nor an LS splitting, we will em- 
ploy spin-averaged masses. For bottomonium, we em- 
ploy the recently measured % mass, 9391(3) MeV, and 
theT mass 9460.3(3) MeV, whose spin average is (M^-j- 
3M r )/4 = 9443(1) MeV. We also have at our disposal 
the P-wave mesons, that are also sometimes believed 
to lie in the regime where perturbative pNRQCD is ap- 
plicable, Xbo w hh mass 9859.4(4) MeV, Xbi with mass 
9892.8(3) MeV, and Xb2 with mass 9912.2(3) MeV, 
that yield a spin average (M Xb0 + 3M Xbl + 5M Xb2 )/9 = 
9899.9(3) MeV. 

Looking towards the 25* bottomonium levels (that, 
although their wave-function starts being immersed in 
the non-perturbative part of the potential, can still to 
some approximation be considered perturbative quarko- 
nium states), we notice that the rjb(2S) mass has not 
been measured yet, so that we take the T(2S) as sole 



Table 2 Experimental data employed to fix the potential pa- 
rameters (adapted from [30]). The error in parenthesis refers 
to the last significant digit. 



{M Vb +3M r )/4 



■3M„ 



'5M X62 )/9 



M 



T(2S) 



(M Vc +3Mj/4,)/4 
M Bc + (3/4) AM h f 



9443(1) MeV 
9899.9(3) MeV 
10023(0.3) MeV 
3067.7(3) MeV 
1.6(4) keV 

6317(8) MeV 



data, without the possibility of spin-averaging. The er- 
ror introduced is expected to be of order 10 MeV or 
less since the hyperfine splitting should be smaller than 
for the ground state (where the difference between the 
T and the spin average is 17 MeV). One shouldn't ex- 
pect this state to be well described in perturbation the- 
ory, but we show results nonetheless for comparison and 
completeness. 

Turning to charmonium, the renowned J /ip state 
has a currently accepted mass of 3096.92(1) MeV, and 
the r) c (lS) a mass of 2980(1) MeV. Their spin averaged 
mass is therefore (M Vc + 3M J/v ,)/4 = 3067.7(3) MeV. 
The NRQCD description of the 2S states in the charmo- 
nium system is generally not accepted, so to obtain one 
more observable at the charmonium scale we turn to 
the S*-wave radiative transition width J/ip — > r/ c j, de- 
cay branch P133 of the J/ijj in the notation of [30]. This 
branching fraction currently yields a width Pjm_hj o7 = 
1.6(4) keV. 

This radiative transition width between the IS" char- 
monium states has also been calculated [31] at NNLO, 
and all we need to do is evaluate it numerically. 

In addition we have at our disposal the mass of the 
B c meson, M Bc = 6277(6) MeV. Here the vector B* is 
not yet known, so that the spin average has to be ob- 
tained by interpolation from other flavor combinations. 
Currently we can lean on the B system for which the 
hyperfine splitting is AM hf = 45.8(3) MeV, the B s 
system where AMhf — 49(2) MeV , and the bottomo- 
nium T — rjb system where AMhf = 69(3) MeV. Our 
interpolated estimate for the B c flavor combination is 
therefore AMhf = 53(3) MeV. Adding this correction 
to the B c mass to yield a spin average Mb c + jAMhj 
gives 6317(8) MeV. 

This parameter set is collected in table 2. 



5.1 Evaluation of Pj/^,_ > ^ c7 

Here we comment on the numerical evaluation of the 
radiative transition between the two IS charmonium 
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states, J/ifj and rj/c. At NNLO we read off [31] 

,3 



r 



_ 16a em 2 k ~t 
= 3 6c M^ X 

A IV1 J/i: 

' x | A a s ((Mj^/2f) 

3 7T 



(35) 



where e c = 2/3 is the charm quark's charge in units of 
the electron charge, and the pole quark mass has been 
eliminated in terms of the physical J/tp mass. Since 
the experimental error in this observable is relatively 
large, it is not sensible to delve into subtleties about 
this elimination. 

The momentum of the photon is given simply as 
fc 7 = (Mj^ - M%J/(2Mj/^) by energy-momentum 
conservation. 

The strong coupling constant is evaluated at two 
different scales in this formula, half the resonance mass 
(corresponding to the scale of the quark mass, or hard 
scale) and the momentum scale of the quark inside the 
resonance pj/^ = r^ 1 (or soft scale). The first one is of 
direct evaluation, but since the characteristic momen- 
tum in charmonium is (p) ~ m r a s = J -^a Sl the second 
occurrence of the coupling constant has a scale that 
depends itself on the coupling constant, a s ( j-^ a^j . 

To obtain a s one then needs to employ Eq. (52) 
recursively. The method is detailed in appendix B.3. 



6 Meson numerical results and parameter fixing 

6.1 Exploratory fits 

In this subsection we gain a feeling for the parameter 
values and explore several alternatives. Table 3 presents 
the result of our fits numbered consecutively, not in- 
cluding the B c mass nor including the 1 jm potential. 

Fits 1, 2 and 3 show the consistency of perturbation 
theory, which is quite reasonable. In these fits the quark 
masses come out consistently in the 1.7 GeV (charm) 
and 5 GeV (bottom) ranges. Here a s (m|), left free, 
varies more significantly upon increasing the order of 
perturbation theory and is largely ensuring that the 
IP bottomonium level is in good agreement with exper- 
iment. In spite of this, the radiative ip — > rj c transition 
width is outside its experimental 2a error band, and 
the splittings in the bottomonium spectrum are signif- 
icantly smaller than experiment. The smallness of this 
splitting pulls the spin-averaged bb(lS) mass to higher 
than physical values. 

The first of these defects is common to cc approaches, 
and within the present scheme it requires certain fine 
tuning. The second problem is related to the fact that 



bottomonium excitations start being sensitive to the 
linear part of the static potential, so that a purely per- 
turbativc quarkonium description is not very precise. 

Fits 4 and 5 show the insensitivity to changing the 
freezing scale for the running coupling constant, that is 
in these fixed to 0.8 GeV, whereas it is 0.5 in all oth- 
ers. Comparing fits 3 and 4 for example, we see that 
the value of this constant is irrelevant for all purposes, 
except a marginal improvement in the bottomonium 
splittings. 

Fit 5 is different in that the coupling constant is not 
allowed to vary, but fixed to the world average at the Z- 
pole. This larger value of the coupling brings the transi- 
tion width to better agreement with experimental data. 
This is due to the NNLO contribution, negative, being 
much enhanced. The charm quark is however pushed to 
the limit of its variation band between 1 and 2 GeV in 
the programme, and increases disagreement with other 
determinations . 

In fits 6 and following we return to a freezing scale 
of 0.5 GeV but change the way to compute the % 2 to 
be minimized. Instead of employing the experimental 



error bands af exp for the quarkonium masses in 



X 



E 



(E th - E ex P)1 



(36) 



we adopt a common error band of 30 MeV for all of 
them. This is in recognition that theory errors for these 
observables are orders of magnitude larger than exper- 
imental errors, and we want to check that the experi- 
mental errors are not weighing the various states unduly 
in the fit. 

In fit 6 we separately fit the two charmonium ob- 
servables and the three bottomonium observables, to 
ascertain the tension between them. This is visible from 
the two different values obtained for the coupling con- 
stant evolved to the Z pole, at the level of 10%. 
In fit 7 we leave the bottomonium 2S and IP exci- 
tations out of the x 2 formula to guarantee that the 
bottom quark mass is fixed to the bottomonium spin- 
averaged ground state. We see that the width rj/^,^ Vc 
pulls the charm quark mass again to the limit of our 
allowed variation band. 

In fit 8 instead we decouple /j/^-^ from the min- 
imization. This immediately relaxes back the charm 
quark mass. 

In fit 9 we use as input the spin-averaged ground 
state masses of bottomonium and charmonium, together 
with the pseudodata a s (m 2 z ) = 0.12, that can be under- 
stood as a fit to the ratio of radiative to total widths 
of the T [21]. At the r pole the running coupling is 
also in agreement with r data, that suggests a s (r) = 
0.330(25) [34]. The radiative width is now below the 
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Table 3 Numerical results and fit parameters. Units are MeV except for the (IS) radiative decay width Fj/^,_,.^ c7 which is 
quoted in keV. The reference experimental data for the various fits can be found on table 2. The quark masses quoted are the 
pole masses and the MS masses at 3 GeV to facilitate comparison. Notes: > at the limit of the allowed fitting range; *fixed 
from the world average value at the Z boson pole; & running coupling separately fit to the charm and bottom data. 



Fit 


Order 




J/l/>->1,.T 




m a {23) 


m b - b (lP) 


a 3 (m A z ) 


m c 


m b 


m c 


m b 


1 


LO 


3068 


2.7 


9477 


9836 


9913 


0.128 


1720 


4980 


820 


4600 


2 


NLO 


3068 


3.0 


9567 


9808 


9908 


0.111 


1740 


5080 


1160 


4780 


3 


NNLO 


3068 


2.6 


9632 


9820 


9901 


0.095 


1700 


5050 


1330 


4900 


4 


NNLO 


3068 


2.6 


9536 


9806 


9911 


0.099 


1720 


5040 


1340 


4860 


5 


NNLO 


3068 


1.9 


9443 


9813 


9986 


0.1184* 


2000* 


5220 


1350 


4970 


6 


NNLO 


3070 


2.4 


9630 


9830 


9910 


0.107/0.097 & 


2000T 


5060 






7 


NNLO 


3067 


2.4 


9443 


9548 


9558 


0.108 


2000+ 


5170 


1530 


4970 


8 


NNLO 


3067 


2.6 


9690 


9824 


9906 


0.095 


1700 


5050 


1368 


4900 


9 


NNLO 


3068 


0.74 


9443 


9685 


9525 


0.12 


2188 


5381 


1540 


5140 


10 


NNLO 


3068 


1.6 


9443 


9688 


9796 


0.111 


1853 


5100 


1340 


4878 



experimental value, showing that with fine tuning of 
the coupling constant it can be brought to the physical 
value. However the IP bottomonium state is now very 
far off, due to the increased coupling, and the quark 
masses are far from other determinations. 

We proceed to fit 10, in which we force the repro- 
duction of the precise experimental number for the IS 
radiative width of 1.6 keV. This happens for a s (m 2 z ) = 
0.111 at NNLO. As expected, this is possible at the ex- 
pense of losing agreement with the IP mass. 

Since the experimental error in the radiative width 
is so much larger than the error in the IP mass mea- 
surement, a best fit will try to compromise by lowering 
the coupling constant in spite of this deteriorating the 
computation of the width. 

Overall it appears that to obtain a perfect value for 
Pj/^i-^Tje requires a little fine tuning, and that in no 
case is it possible to obtain an excellent fit to all five 
quantities simultaneously. 



6.2 Extended fits 

In this subsection we include the B c mass and explore 
in addition the effect of the 1/m potential. In all cases 
the coupling constant is evolved to the Z pole at NNLO 
(this should bc for broad comparison and not taken as a 
detailed prediction since higher orders of perturbation 
theory should then be used). 

Comparing tables 4 and 5 we see that the effect of 
the 1/m potential is modest, the fit preferring a slightly 
lower coupling constant, and the B c mass being better 
adjusted. 

Tables 6 and 7 then show the same calculation but 
in the PS scheme. An interesting feature in these com- 
putations is seen in the last three rows of table 7. When 
the 1/m correction is included, the PS scheme seems to 



Table 4 Further fits in pole scheme. Masses are in MeV , the 
radiative decay width J/tp — > r/c'y in keV. The infrared freez- 
ing of the running coupling constant occurs at 0.6 GeV. The 
recoil 1/m potential is not included, only the static potential. 





LO 


NLO 


NNLO 


Expt. 


Fit number 


1 


2 


3 






3068 


3068 


3068 


3068 


1 J / il>—>-r} c ~Y 


2.7 


3.1 


2.1 


1.6(4) 


m bb (lS) 


9458 


9447 


9480 


9443 


m bb (2S) 


9820 


9777 


9769 


10023 


m bb (lP) 


9899 


9900 


9897 


9900 


m bc (lS) 


5922 


6158 


6158 


6317 


m c 


1720 


1770 


1850 




m b 


4970 


5040 


5120 




a a 


0.128 


0.116 


0.111 





Table 5 As in table 4 but with the recoil 1/m potential 
included. 





LO 


NLO 


NNLO 


Expt. 


Fit number 


4 


5 


6 






3068 


3068 


3068 


3068 


1 J / il>—>-r} c ~Y 


2.7 


3.0 


2.5 


1.6(4) 


m bb (lS) 


9443 


9445 


9488 


9443 


m bb (2S) 


9775 


9788 


9772 


10023 


m bb (lP) 


9900 


9900 


9896 


9900 


m bc (lS) 


6330 


6210 


6364 


6317 


m c 


1740 


1740 


1820 




m b 


5040 


5000 


5110 






0.102 


0.106 


0.105 





be rather stable in going from LO to NLO to NNLO, 
as the quark masses barely change. 

Comparing tables 8 and 9 we see that, with a lower 
infrared cutoff, the PS scheme does somewhat better in 
terms of convergence and agreement with data. 

Fixing the running coupling constant to either the 
world average or recent r data, as in tables 10 and 11 
leads to slightly improved agreement with experiment 
in the NNLO computation of the radiative transition 
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Table 6 As in table 4 but in the PS scheme with the poten- 
tial totally cutoff in the infrared at 0.6 GeV. 





LO 


NLO 


NNLO 


Expt. 


Fit number 


7 


8 


9 






3068 


3068 


3068 


3068 


P.J / 1/j — ¥7) c "i 


2.7 


3.0 


2.5 


1.6(4) 


m b - b (lS) 


9462 


9443 


9443 


9443 


m b - b {2S) 


9818 


9792 


9789 


10023 


m bb (lP) 


9898 


9900 


9900 


9900 


m bc (lS) 


5866 


6264 


6270 


6317 


m c 


1740 


1930 


1710 




m b 


4970 


5230 


5020 






0.133 


0.112 


0.106 





Table 9 Same as table 8 but in the PS scheme, with the 
potential completely cutoff in the infrared at 0.4 GeV. 





LO 


NLO 


NNLO 


Expt. 


Fit number 


7a 


8a 


9a 






3068 


3068 


3068 


3068 


Pj/ip — ► T] c 7 


2.7 


3.0 


2.5 


1.6(4) 


m bb {lS) 


9614 


9443 


9443 


9443 


m bb {2S) 


9836 


9791 


9788 


10023 


m bb (lP) 


9884 


9900 


9900 


9900 


m bc (lS) 


6218 


6265 


6269 


6317 


m c 


1631 


1718 


1748 




m b 


4964 


5029 


5055 






0.115 


0.112 


0.105 





Table 7 As in table 4 but in the PS scheme with the poten- 
tial totally cutoff in the infrared at 0.6 GeV, and with the 
1/m recoil potential included . 





LO 


NLO 


NNLO 


Expt. 


Fit number 


7 


8 


9 






3068 


3068 


3068 


3068 


Pj/ i/j— 7-77^7 


2.6 


3.0 


2.5 


1.6(4) 


m bb (lS) 


9443 


9444 


9443 


9443 


m bb {2S) 


9791 


9792 


9789 


10023 


m bb {lP) 


9900 


9900 


9900 


9900 


m bc (lS) 


6270 


6264 


6270 


6317 


m c 


1690 


1700 


1710 




m b 


5000 


5020 


5010 




a s 


0.106 


0.112 


0.106 





Table 8 Further fits. Masses are in MeV, the radiative decay 
width J/%p — > r? c 7 in keV . The infrared freezing of the running 
coupling constant occurs at 0.4 GeV . 





LO 


NLO 


NNLO 


Expt. 


Fit number 


la 


2a 


3a 






3068 


3068 


3068 


3068 


Pj/ Tp^?) c "{ 


2.7 


3.1 


2.5 


1.6(4) 


m bb (lS) 


9542 


9481 


9577 


9443 


m bb (2S) 


9829 


9800 


9804 


10023 


m bb (lP) 


9891 


9897 


9888 


9900 


m bc (lS) 


6112 


6118 


6277 


6317 


m c 


1667 


1828 


1761 




m b 


4963 


5117 


5084 




a s 


0.120 


0.121 


0.104 





Table 10 Further fits fixing the coupling constant at the Z- 
pole at cts{mz) — 0.1204 as determined by recent r data. 
Masses are in MeV , the radiative decay width J/tj; —> r/cf in 
keV. The infrared freezing of the running coupling constant 
occurs at 0.4 GeV. 





LO 


NLO 


NNLO 


Expt. 


Fit number 


lb 


2b 


3b 






3068 


3068 


3068 


3068 


p 

L J 1 'ip—*ri c ~f 


2.7 


3.0 


0.65 


1.6(4) 


m bb (lS) 


9444 


9443 


9444 


9443 


m bb (2S) 


9669 


9721 


9586 


10023 


m hb {lP) 


9718 


9812 


9654 


9900 


m bc (lS) 


6145 


6195 


6149 


6317 


m c 


1629 


1728 


1889 




m b 


4872 


5000 


5116 




Q-s 








fixed at 0.1204 



Table 11 Same as table 10 but in the PS scheme, with the 
coupling constant cutoff at 0.4 GeV. 



Fit number 


LO 
7b 


NLO 

8b 


NNLO 
9b 


Expt. 


nice 


3068 


3068 


3068 


3068 


Pj/ ip— *-f? c 7 


2.7 


3.0 


0.65 


1.6(4) 


m bb (lS) 


9444 


9443 


9443 


9443 


m bb (2S) 


9646 


9791 


9885 


10023 


m bb {lP) 


9690 


9900 


10034 


9900 


m bc (lS) 


6155 


6264 


6248 


6317 


m c 
m b 
a s 


1619 
4855 


1718 
5029 


1918 
5243 


fixed at 0.1204 



width of the J/?p, but in exchange the IP bottomonium 
mass is completely off. 

The difference between table 10 and entry number 
9 of table 3 is that here the global fit strategy was used 
while there only the IS" masses were used to constrain 
the quark masses, and that the freezing of the coupling 
occurs at a slightly different momentum (0.4 versus 0.46 
GeV). 

The charm quark mass has recently [33] been reob- 
tained from a lattice computation of the ground state 
charmed and charmonium mesons (D, J/ip, i] c ), with 



Nf = 2. In the MS scheme they obtain 

mf^(2 GeV) = 1.14(4) GeV (37) 

that translates into a mass at the charm scale of 

m c (m c ) = 1.28(4) GeV . 

In comparing the quark mass between various schemes 
and scales [35], the collaboration quotes an error less 
than one standard deviation as a result of using Nf = 2 
instead of Nf = 4. 
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Table 13 Parameters employed by Yu Jia in an early varia- 
tional computation of the triply heavy baryon spectrum with 
a strong Coulomb potential. The coupling constant is given 
at the charm scale (employed in the ccc computation) and 
the bottom soft scale (employed for all others). 

q s (0.9 GeV) (159 
a s (1.2 GeV) 0.43 

mc = Mm.(i + 2^ ~ 1.668 GeV 
m b = ^(l + 2^) - 4.924 GeV 



Translating our pole-scheme masses in the various 
tables into MS masses m c (m c ) consistently yields re- 
sults of order 1.3 GeV, in agreement with the lattice 
determination. The PS scheme on the other hand gives 
a somewhat larger mass of about 1.4 GeV in the MS 
scheme. 

To conclude this section let us quote an NLO com- 
putation with the BLM scheme, with best fit a s = 
0.113, m c = 1750 MeV and m b = 5030 MeV. The cc 
radiative width is 3.0 keV. The bb (1S,2S,1P) masses are 
9393, 9775 and 9904 MeV respectively. The B c mass 
comes out to be 6180 MeV . Thus, there is no particular 
advantage in using this scheme over the PS or the pole 
ones. 



7 Prior computations of the baryon masses 

In this section we compile existing computations of the 
various ground state triply heavy baryons. Since we 
work in the static limit, we will not resolve the hyper- 
fine spin splitting between spin 1/2 and spin 3/2. Thus, 
in table 12 we quote the spin average (M 1 / 2 + 2M 3 / 2 )/3 
for the bbc and ccb wave-functions, that can appear in 
both spin combinations in the ground state. We further 
plot some of these computations in figure 4. 

Most computations in the literature are consistent 
with a triply charmed baryon of around 4800 MeV, 
and ours will not bc different. The bbb baryon is pre- 
ferred by most approaches in the range of 14400 MeV. 
A salient exception is the sum-rule computation of [3G] 
that seems to be significantly lower. 

Closest in spirit to our approach is the Coulombic 
calculation of Jia [37], that could be considered a Lead- 
ing Order pNRQCD computation of the static poten- 
tial. Indeed, this author employs Eq. (21), with param- 
eters specified in table 13. The error bar quoted in that 
work corresponds to the author's estimate of higher or- 
ders in perturbation theory. We will explore the system- 
atics of this computation, extending it in several ways. 
First, we will work to two higher orders in perturba- 
tion theory with the potentials now available. Thus, 



we will ascertain that this error was underestimated. 
Second, we will quantify the error implicit in the one- 
wavefunction variational approximation (that Jia also 
uses) by showing explicit computations in very similar 
atomic systems for which the experimental data is avail- 
able. And third, we will incorporate a running coupling 
constant at all steps, and handle the attending infrared 
systematic uncertainties by comparing different meth- 
ods. The outcome of our work will thus be a much more 
detailed understanding of triply heavy baryons in the 
context of pNRQCD. 

It is also worth remarking that in [38] the difference 
between the triply heavy baryon mass with a zi-like 
two-body potential and an T type potential has been 
reported in a variational model computation. We have 
plotted their results in figure 5. As is easily seen, those 
authors find that the T configuration is slightly heavier, 
but in any case the difference is only of order 20-40 
MeV. 

This observation is of interest for light-quark baryons, 
since they can be more conveniently treated by covari- 
ant means, but the (much simplified) Faddeev equations 
require a vanishing pure three-body kernel. That there 
is not much difference between A and T shape config- 
urations is important information for establishing that 
the Faddeev equations are approximately valid, at least 
in the heavy quark limit when the soft scale is pertur- 
bative 4 . 

In our perturbative treatment for heavy baryons we 
will compare the mass computation with and without 
the pure three-body force that appears at NNLO in 
perturbation theory, finding that this difference is also 
small, and thus further reinforcing the conclusion of 
Flynn et al.. 

8 Novel computation of triply heavy baryon 
masses 

We treat the 3-body problem in a similar manner to [53, 
54] variationally by employing a simple wave-function 
ansatz and computing the expectation value of the pN- 
RQCD Hamiltonian. The Raylcigh-Ritz variational prin- 
ciple guarantees that the outcome is an upper bound of 
the true ground state energy in the given channel, 

{^a p a x \Hp NR QCD\^ apa> ) , . 

Ti ui \ ' ( > 

The two parameters a p , a\ are then varied to find 
the best possible upper bound on energy for the given 

4 Should the soft-scale be non-perturbative, one could ob- 
tain some information from other published work[51,52], but 
incorporating this in a three-body computation is currently 
beyond our scope 
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Table 12 Computations of triply heavy baryon masses. Many of the entries were already compared in [38] 
more information available in the literature. 



but we have added 



Method 


Ref. 


M ccc 


M cc f, 


M bbc 


Mbbb 


Variational Coulomb 


[37] 


4760(60) 


7980(70) 


11190(80) 


14370(80) 


Variational Cornell 


[38] 


4799 


8037 


11235 


14398 


Faddeev 


[39] 


4799 


8019 


11217 


14398 


Bag model 


[40] 


4790 


8030 


11200 


14300 


Quark counting rules 


[43] 


4925(90) 


8200(90) 


11480(120) 


14760(180) 


Const, quark model 


[44] 


4965 


8258 


11548 


14834 


Const, quark model 


[45] 


4632 








Relat. quark model 


[46] 


4803 


8023 


11285 


14569 


Instanton quark model 


[47] 


4773 








Hypercentral model 


[48] 


4736 


8096 


11381 


14451 


Sum rules 


[36] 


4670(150) 


7443(150) 


10460(110) 


13280(100) 


Lattice 


[49] 


4780 






14371(12) 


Regge estimate 


[50] 


4819(7) 









■4 ► # 



# Hasenfratz'80 (Bag model) 

■ Bjorken 

y Vijande'04 (NRQM) 

▲ Roberts'08 (NRQM) 

-4 Martynenko'08 (RTQM) 

+ Flynn' 1 1 (NRQM) 

$ Silvestre-Brac'96 (Faddeev) 

x Guo'08 (Regge) 

© Zhang'()9 (Sum rules) 

0—0 Jia'06 (Coulomb, variational) 

+ McNeile' 1 1 (lattice) 



15000 c— 
14900 i- 
14800 =- 
14700 =- 
14600 j- 
14500 E- 
14400 j- 
14300 =- 
14200 =- 
14100 j- 
14000 '- 
13900 z - 
13800 =- 
13700 =- 
13600 L 
13500 %■ 
13400 =- 
13300 %■ 
13200 =- 
13100 =- 
13000 



► * 



• Hasenfratz'80 (Bag model) 

■ Bjorken 

▲ Roberts'08 (NRQM) 

■4 Martynenko'08 (RTQM) 

k. Flynn' 1 1 (NRQM) 

jK Silvestre-Brac'96 (Faddeev) 

© Zhang'09 (Sum rules) 

3— H Jia'06 (Coulomb, variational) 

+ Meinel' 10 (lattice) 



Fig. 4 Scatter of several existing computations for the Q ccc and fl^bb masses respectively. See table 12 for references. 




V 



I 

# 



• 


ccc 


■ 


ccb 


▲ 


bbc 


♦ 


bbb 



I 

0.1 



Fig. 5 Mass difference between triply heavy baryons computed with the A type potential (two-body interactions alone) and 
T-type potential (three-body interactions alone) reported by Flynn, Hernandez and Nieves in [38]. 



ansatz. These two parameters are associated to two 
momentum-space Jacobi coordinate vectors, the third 
independent vector being fixed by the center of mass 



condition (hadron at rest) 
k\ - k 2 



hp — 



kx = V 2^ + k ^ 
k 3 = —ki - k 2 ■ 



V2 
3 



(39) 
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Table 14 Ground state triply heavy baryon masses in the 
PS scheme, with infrared cutoff A = 0.4 GeV, for various 
orders of perturbation theory. All masses in MeV . 



Order 


ccc 


ccb 


bbc 


bbb 


LO 


4895 


8235 


11535 


14770 


NLO 


5160 


8480 


11750 


14970 


NNLO 


5250 


8560 


11805 


15040 



Table 15 Ground state triply heavy baryon masses in the 
PS scheme, with infrared cutoff A = 0.6 GeV, for various 
orders of perturbation theory. All masses in MeV. 



Order 


ccc 


ccb 


bbc 


bbb 


LO 


5240 


8500 


11640 


14750 


NLO 


5810 


9170 


12460 


15670 


NNLO 


5150 


8690 


12100 


15500 



We choose as ansatz 

i>(k p , k x ) apax = Y m (k p )Y o(kx)e- k > /a ?- k > /a * , (40) 

which gives reasonable results (we have also checked 
other forms such as a rational function). The error in- 
curred in this variational approximation is estimated 
below in subsection 9. The wave-function in Eq. (40) 
is symmetrized as needed by invoking it in the com- 
puter programme with different arguments, to guaran- 
tee symmetry under exchange of any two equal quarks. 
The color singlet wave- function e^fc /a/3 that is implicit 
in the calculation (and has already been used in the 
computation of the color factors of the various poten- 
tials) is then responsible for the antisymmetry expected 
under Fermion exchange. 

In practice we compute the Hamiltonian's expecta- 
tion value for the three-body problem in momentum 
space. 

8.1 Results in the PS scheme 

With the three parameters a s , m c , in the PS scheme 
fit to the meson spectrum, the only sensitivity left to 
explore is that of the infrared cutoff scale. In tables 14 
and 15 we present the outcome of the three-body com- 
putation in the PS scheme. 

The mass values obtained are significantly higher 
than in other approaches. As will be seen in the next 
section, this is a feature of the PS scheme, that misses 
quite some of the binding, as opposed to the pole scheme. 
This feeling is reinforced by the observation that the 
computation with the lower 0.4 GeV cutoff does much 
better, both in terms of binding and convergence. Par- 
ticularly bad is the computation with an infrared cutoff 
at 0.6 GeV at NLO, that yields an unbelievably high 
mass. 



Table 16 Ground state triply heavy baryon masses in the 
Pole scheme, with infrared freezing point A = 0.4 GeV , for 
various orders of perturbation theory. All masses in MeV . 



Order 


ccc 


ccb 


bbc 


bbb 


LO 


4708 


7975 


11180 


14386 


NLO 


4900 


8140 


10890 


14500 


NNLO 


4865 


8150 


11400 


14683 



Table 17 The difference between the PS scheme and the pole 
scheme seems to persist at masses twice and thrice as big as 
the bottom quark, with the PS scheme underestimating the 
binding energy. Although for asymptotically large masses we 
believe that this difference should ameliorate, we do not see 
it presently. 



Scheme 


Quark mass (GeV) 


Baryon mass (GeV) 


Pole 


5.08 


14.68 


PS 


5.06 


15.04 


Pole 


10 


29.22 


PS 


10 


29.67 


Pole 


15 


43.95 


PS 


15 


44.46 



The results in this and the next subsection satisfy 
Nussinov's inequalities [41]. The first, 

M Qbbc < 2Ma aei - M Qcac (41) 

is a consequence of heavier systems being more bound 
than lighter systems (as discussed at the end of subsec- 
tion 2.1). The second inequality, satisfied by a sizeable 
amount, reads 

Mr 

M fifc >y+% (42) 

and means that mesons are more tightly bound than 
baryons. They are well satisfied when the three-body 
computation is compared to the corresponding two- 
body computation under the same scheme and condi- 
tions employed for parameter fitting. 



8.2 Results in the pole scheme 

Tables 16 and 18 present our results in the pole scheme 
with couplng freezing at 0.6 and 0.4 GeV respectively. 

Comparing tables 16 and 14 we see that the PS 
scheme, with its drastic infrared cutoff to avoid renor- 
malons, is underestimating the binding energy by a 
large amount of order 300 MeV. To check whether this 
is ameliorated for yet heavier quarks we have ran also 
with 771q = 10 GeV, a quark twice as heavy as the bot- 
tom, and with niQ — 15 GeV. The results are shown in 
table 17. 
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Table 18 Ground state triply heavy baryon masses in the 
Pole scheme, with infrared freezing point A = 0.6 GeV, for 
various orders of perturbation theory. All masses in MeV. 



Order 


ccc 


ccb 


bbc 


bbb 


LO 


4750 


7950 


11100 


14200 


NLO 


5050 


8290 


11470 


14630 


NNLO 


4970 


8200 


11340 


14570 



H-H Jia'06 (Coulomb, variational) 
• LO pNRQCD 
■ NLO pNRQCD 
Y NNLO pNRQCD (A pot. only) 



5800 
5700 
5600 
55O0 
54O0 
5300 
5200 
51(1(1 
5000 
4900 
4800 
4700 
4600 
450(1 



Fig. 6 Q ccc computations in Pole scheme (solid) and PS 
scheme (shaded, red online) to LO, NLO and NNLO. The 
leftmost point (green online) is the original Coulomb evalu- 
ation by Yu Jia with his quoted error band estimating the 
NLO effect. Infrared saturation (in the pole scheme) or cut- 
off (in the PS scheme) set at 0.6 GeV . The asymmetric error 
band is our extrapolation of the missing binding energy due 
to the variational wave-function. 



For the 10 GeV quark we obtain a mass of 29.22 
GeV in the pole and 29.67 GeV in the PS scheme. Al- 
though the difference is now a smaller percentage of the 
total mass, it is still very significant in absolute terms. 
Therefore we do not expect that a small refitting of pa- 
rameters, such as quark masses or coupling constant, 
will be able to eliminate it within the present setup. 

After all numbers have been considered, we deem 
that the computation that has the best balance between 
convergence of perturbation theory and capture of the 
infrared physics is that in table 16. To obtain the best 
estimate of the physical baryon mass, the results com- 
puted there have to be extrapolated by increasing the 
binding energy by 25% to compensate for the varia- 
tional approximation (see section 9). 

Figure 6 compares the results in the PS and pole 
schemes with an infrared saturation scale of 0.6 GeV , 
to the three orders of perturbation theory available. 

Several conclusions follow from the figure. First, it is 
plain that Jia's calculation is in the right ballpark, but 
underestimates the corrections due to higher orders of 
perturbation theory (note that our coupling constant is 
on the low side of the world average, such that a scheme 
that will reduce these corrections is hard to imagine). 
In addition, one sees that as already mentioned, the PS 



3— H Jia'06 (Coulomb, variational) 
• LO pNRQCD 
■ NLO pNRQCD 
y NNLO pNRQCD (A pot. only) 



Fig. 7 Same as in figure 6 but for the 1R saturation/cutoff 
at a lower scale of 0.4 GeV. 



scheme underestimates the binding. Finally, and taking 
into account the variational error bar (any such calcu- 
lation underestimates the binding), the prediction for 
the ccc mass should be about 4800 MeV. We later will 
correct this figure up when accounting for the V m -i 
potential in subsection 8.4. 

Comparing with the results plotted in figure 7, we 
see that the region between 0.6 and 0.4 GeV still con- 
tributes at least an additional 100 MeV of binding. 

Figure 8 then plots the predictions for the other 
(spin-averaged) triply heavy baryons (in the pole scheme) 
and gives a panoramic of other results in the literature. 

Figure 9 shows the size of the binding energy from 
the three body variational calculation by comparing it 
with three times the pole mass (given that the potential 
we use is extracted from perturbation theory, this acts 
as a dissociation threshold, that should not be present 
in a lattice or a Cornell model computation, for exam- 
ple). It is plain that, although separately the pole mass 
and the binding energy do not converge well, there is 
a cancellation between them that helps the behavior of 
the baryon mass in perturbation theory. 



8.3 Effect of the three-body force 

Next we address the difference between a computation 
employing the intrinsic three-body force, and a compu- 
tation with only the two body force. If we set the scale 
in Eq. (26) according to the "hard scale" prescription 

a 3 s -> a s (m c ) 3 

we obtain a very small pure three-body contribution of 
order 1 MeV. If instead, in view of the typical momen- 
tum transfer through the three gluons, we choose the 
more sensible "soft scale" prescription 



a 3 s -> a s (q 2 ) a s (q 3 ) a s (y/q 2 q$ ) 
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► # 



• Hasenfratz'80 (Bag model) 
■ Bjorken 

* pNRQCD (LO.NLO.NNLO) 
▲ Roberts'08 (NRQM) 

4 Martynenko'08 (RTQM) 

Flynn' 1 1 (NRQM) 
$ Silvestre-Brac'96 (Faddeev) 
E3— H Jia'06 (Coulomb, variational) 
+ Meine!' 10 (lattice) 



Hasenfratz'80 (Bag model) 
Bjorken 

pNRQCD (LO,NLO,NNLO) 
Roberts'08 (NRQM) 
Martynenko'08 (RTQM) 
Flynn' 11 (NRQM) 



H— H Jia'06 (Coulomb, variational) 



• Hasenfratz'80 (Bag model) 
■ Bjorken 

* pNRQCD (LO.NLO.NNLO) 
▲ Roberts'08 (NRQM) 

4 Martynenko'08 (RTQM) 

^ Flynn' 1 1 (NRQM) 

EH— H Jia'06 (Coulomb, variational) 



Fig. 8 Predictions for the mass of the Qbbb, and spin aver- 
aged ccb, bbc combining results analogous to those of figure 6 
and including other computations as in figure 4. 
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15200 
15100 
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14900 
14800 
14700 
14600 
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14400 
14300 
14200 
14100 
14000 

Fig. 9 Binding energy of the triple-6 system. Lines: three 
times the pole mass (from bottom to top LO, NLO, NNLO). 
Symbols: mass of the Qbbb- The graph shows how the res- 
onance mass is better behaved in perturbation theory than 
either of the pole mass or the binding energy. 



. . 3M b LO 

3M b NLO 

» pNRQCD (LO,NLO,NNLO) 
3M, NNLO 



• 


ccc 


■ 


ccb 


♦ 


bbc 


▲ 


bbb 



Fig. 10 The effect of adding the three body force to the 
NNLO potential is to increase the binding in the amounts vis- 
ible. The effect is larger for equal flavor objects (more tightly 
bound) than for mixed flavor, and does not depend much on 
quark masses. 



.4 Effect of the 1/m potential 



the effect of the (perturbative) three-body force is of 
order 20-40 MeV, in broad agreement with the related 
(though not equivalent) estimates of Flynn et al. The 
result of the computation of the three-body potential 
for the different flavor combinations is depicted in fig- 
ure 10. The effect we find is of size 17 MeV for ccc, 25 
MeV for ccb, 39 MeV for bbc and 20 MeV for bbb, 
with an error of about 5 M eV or less. 

The immediate conclusion is that in the heavy quark 
limit, intrinsic three body forces (defined as those van- 
ishing when one quark is put far away from the other 
two) are small in ground state baryons. 



Thus far our three-quark results have been based on the 
purely static potential. In this section we lift this ap- 
proximation and study, at NLO, the effect of adding a 
V m -i contribution. This recoil correction has not been 
worked out in detail in the literature, so we abstain 
from attempting an NNLO evaluation. But if we turn 
to the simplest convention of [7], the NLO V m ~i is en- 
tirely given by the non-Abelian diagram with a three 
gluon vertex, whose equivalent for baryons is sketched 
in figure 11. 

Because the interaction is two-body, the potential in 
Eq. (13) can immediately be adopted for baryons, with 
appropriate kinematics and excepting a color factor. 
The latter can be worked out easily by noting that the 
diagram is a one-loop radiative correction to the quark- 
gluon vertex on the quark at the very top of figure 11. 
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Fig. 11 Non-Abelian diagram that produces a recoil V m -i 
potential in triply heavy baryons at NLO. 



This non-Abelian vertex correction is easily seen to in- 
crease the usual Gell-Mann color matrix at the ver- 
tex [42] 



300 



250 h 

> 

§ 

H 200 
< 



150 



• 


ccc 




■ 


ccb 




♦ 


bbc 






bbb 





100 

Fig. 12 Effect of V m -i at NLO from Eq. (13) on the ground- 
state baryon spectrum. Plotted is Mi — Mo, the mass dif- 
ference including the 1/m potential or employing the static 
potential alone. The mass of the Q ccc can be raised by more 
than 150 MeV due to recoil corrections. 



rjia y c rpg 

2 

Thus the ratio of color factor in baryons over color fac- 
tor in mesons is the same for the static potential as 
for the recoil correction, to NLO and in the convention 
of [71 



a 



3.m 



c. 



3,0 



Co 



C2,0 



(43) 



and in practice it is sufficient to divide Eq. (13) by a 
factor 2 to obtain each of the three possible two-body 
interactions in the baryon system. 

Next we will show the difference in baryon mass 
with and without the V m -i potential. To properly nor- 
malize the pole mass and coupling we first recompute 
the meson spectrum in section 6 and ensure a best fit 
shown in table 19. 

To be consistent with the given order in perturba- 
tion theory, the coupling constant runs only at NLO. 
The quark mass takes an almost identical renormaliza- 
tion of —56 MeV (for charm) or —61 MeV (for bottom) 
upon including V\i m that carries over to the baryon 
computation and is accounted for in addition to the 
recoil interaction there. 

The difference between computing baryon masses 
with the recoil potential or without it at NLO is de- 
picted in figure 12. As can be seen, the effect increases 
softly from ccc (194(3) MeV) to bbb (297(3) MeV) as 
discussed earlier around Eq. (16). 



9 Error estimates 

All integrals in the computation of the three-body Hamil- 
tonian matrix element are evaluated by Monte-Carlo 





Para-Helium 



Ortho-Helium 



Dihydrogen cation 



Fig. 13 Three-body systems in atomic physics that we use 
to test the variational method and computer programme. 
From left to right, para-Helium (one a particle and two elec- 
trons with spin antialigned), ortho- Helium (the two electron 
spins are now aligned) and the Hydrogen cation (two protons 
loosely bound by one electron alone). 



methods and we allow an error of 10 MeV in their com- 
putation, except in our three-body force or recoil force 
computations. In those we have demanded an error in 
the 1-5 MeV range given that we have to subtract two 
masses. This numerical uncertainty will be negligible in 
the final error balance. 

To estimate the variational errors we turn to some 
simple systems in atomic and molecular physics that 
can be addressed with the same techniques, providing 
in addition a check of the computer programmes. We 
take three-body systems made of one electron and two 
protons (the dihydrogen cation ), and one a-particle 
binding two electrons with parallel or antiparallel spins 
(ortho and para-Helium respectively). These are de- 
picted in figure 13. 

Although we are not considering spin interactions, 
the distinction between ortho and para-Helium is also 
important as it checks our wave-function symmctriza- 
tion procedure. 

We give the matrix elements in terms of reduced 
momenta k = k/{ra e a em ). They read, for atomic He- 
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Table 19 Meson observables at NLO in the pole scheme, with coupling constant freezing only at the very low 0.4 GeV scale. 
All (spin averaged) masses in MeV. 



Potential 


cc(lS) 


66(15) 


66(1P) 


B c 


ti s (mz) 


m c 


m b 


V m -i 


3068 


9443 


9914 


6200 


0.107 


1767 


5032 


Static 


3068 


9443 


9857 


6104 


0.121 


1823 


5093 



Para-He -79.0 eV (calc. -59eV) 
Ortho-He -59.2 eV (calc. -48eV) 

H, + -15.4 eV (calc. -12.5eV) ♦ 



Or 
-10- 
-20- 
-30- 
|-40 L 

-60- • ■ 

-70- 

-80 - • 

-90 1 

Fig. 14 Variational estimate of the various binding energies 
in atomic three-body problems, together with the experimen- 
tal values. 



Hum 

Itt\ 2 f d 3 kl rf3fc 2 ,»,, , , 

{H)^ = m e a em J (27r)3 (27r)3 ^ {h,k 2 )x 



(44) 



AirtPq 
(2tt)V 



N)(ki + q,k 2 - q) - 2tp(ki + q, k 2 ) 
-2ip(k 1 ,k 2 +q) 



and for the diHydrogen cation 
d 3 /ci d 3 k 2 



(H)^ = m e a 



(2tt) 3 (2tt) 3 



1 



2^ 1 + (kl + kl)^Ji>(h,k 2 )+ (45) 

(27r)V ( ~ + <7 ' &2 ~ ^ ~ + 9 ' ^ 

+ V(*l.*2 + «)) 



where the 4ir/q 2 Coulomb potential is clearly recogniz- 
able, and makes the expression very much alike to our 
heavy-baryon computation at LO. 

The outcome of these atomic computations is then 
plotted in hgure 14. 

From this exercise we estimate the precision of the 
simple one-wavefunction, two-parameter variational es- 
timate, to bc about 25% in error in the computation 



of the binding energy. This error can eventually be re- 
duced to zero by a systematic shell by shell diagonal- 
ization with a large basis, and we arc planning such 
undertaking, but it exceeds the purposes of the present 
work. 

For now we just note that since the Rayleigh-Ritz 
variational principle guarantees that the computed state 
is less bound than the physical state, we can reduce the 
error by the device of increasing the binding energy 
that estimated 25% in the final quoted estimate. This 
extrapolation is assisted by our atomic physics compu- 
tation, since once the respective scales are pulled out 
of the matrix elements for heavy baryons or for atomic 
Helium, they are very similar. 

Returning next to heavy baryons, we observe that 
perturbation theory seems to be converging reasonably 
well, and that, while the jump from LO to NLO is ap- 
preciable, the difference between NLO and NNLO is 
substantially smaller and of order 100 MeV at most. 
This reasoning applies also to higher order recoil cor- 
rections. 

We also incur in a small inaccuracy of order 5 — 10 
MeV in the computation of the ccb and bbc mixed- 
flavor mesons in employing Nf = 3 in the appropriate 
coupling constant, instead of varying the screening Nf 
with the scale of the various interactions, that may re- 
quire one of Nf = 3 or Nf = 4. This is in order to sim- 
plify and speed the execution of computer programmes. 
In neglecting the charm sea in these baryons we are in 
line with many other modern computations. Nf = 4 is 
correctly set for bbb baryons. 

The error is therefore dominated by the treatment 
of the infrared. The Potential Subtracted Scheme offers 
perhaps slightly improved convergence in perturbation 
theory, but since the potential is completely truncated 
at a low scale, it underestimates the binding energy 
systematically, and therefore overestimates the mass by 
a significant amount (several hundred MeV). 

This notorious effect should be absent for infinitely 
heavy quarks, where mv 3> Aqcd, but since for physi- 
cal charm and bottom quarks the scale separation is not 
very clean, we see that imposing an infrared cutoff as 
the PS scheme demands affects the computed binding 
energy. 

Turning to the pole scheme with saturated running 
constant in the infrared, we see that the binding energy 
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is more in line with other approaches, and the conver- 
gence of perturbation theory is still acceptable. How- 
ever the difference between the two approaches advises 
us to assign an error of order 200 MeV to the result. 

That this infrared sensitivity is caused by the finite- 
ness of the quark mass can be exposed by the following 
reasoning. In the case of the Hydrogen atom, the bind- 
ing caused by momenta smaller than the inverse Bohr 
radius 1/ao is about 13% of the total —13.6 eV. In 
QCD the running coupling constant continues growing 
below Oq 1 which is of order 500 MeV, causing a corre- 
spondingly higher error. For quarks of much larger mass 
however, Oq 1 is higher and the low momentum tail of 
the wave-function has little overlap with the infrared 
region where ^4qcd influences the result. 



10 Prospects for experimental detection 

James Bjorkcn [43] proposed several decay chains ac- 
cessible to experiment that would allow the reconstruc- 
tion of the ft C cc- More recently, Chen and Wu [55] have 
estimated that 10/6 -1 of LHC integrated luminosity 
recorded by a detector would contain 10 4 to 10 5 triply 
charmed baryons. They further propose looking for the 
particular decay chain 

i? ccc —> ftsssiririr . 

We agree that this channel provides in principle a very 
clean signature, since all four final state particles are 
charged and can be identified by dE/dx energy deposi- 
tion on a tracking chamber. However, the tremendous 
combinatorial background that the LHC experiments 
have to contend with, often with hundreds of pions in 
a single event, make the search extremely difficult. 

The reconstruction could be performed also at other 
experiments, such as perhaps COMPASS [57] and cer- 
tainly the B factories, operating at less energy and thus 
with less multi-particle production. These, particularly 
COMPASS, are however limited in statistics and pro- 
duce less triple charm events. Therefore, it would be 
an advantage to combine the spectra in all Cabibbo- 
allowed four charged particle channels f2 sss innr, SKtttt, 
UK Kir, and pKKK. 

Still, an alternative route to complete reconstruc- 
tion of the fi CC c would be to at least measure its mass 
in a recoil spectrum. Since charm is produced by the 
strong interactions in cc pairs, the f2 ccc needs to recoil 
against three charm antiquarks, most often in the form 
of three D mesons, for example in the reaction 



e~e+ -> Q ccc p DDD (46) 



that, depending on the energy, can be accompanied by 
any number of pions. 

The interesting feature to exploit is that the C cc 
is the lightest state that can carry three charm quarks. 
Thus, an alternative strategy to reconstructing a decay 
chain, any of which will have small branching fractions, 
would be to search for the recoiling triple (anti)charm 
system and recoil antiproton. 

Since no triple charm spectrum has ever been pub- 
lished to our knowledge, searching for three D mesons is 
an interesting undertaking in itself. This can be accom- 
plished, for example, by a lepton trigger that tags one 
of the charm mesons (only 50% inefficient while sup- 
pressing much background), followed by reconstruction 
of the other two, or by a pure hadronic trigger in which 
all three are fully reconstructed. 

The further identification of a recoil antiproton in 
a small subset of the events immediately provides an 
upper bound on the i? ccc mass, simply by the missing 
energy technique against the recoiling system, even if 
the fl ccc itself was not fully reconstructed. 

Baranov and Slad [56] estimated that the produc- 
tion cross section for fi ccc DDD at the Z-polc in e~e + 
collisions is of order 0.04/6 -1 , which is too small to have 
been usable at LEP or SLC. However, we should take 
into account that the i flux factor allows for a larger 
cross-section, up to perhaps 3fb, in the 10 GeV region 
where the B factories operate. Such cross-section is not 
unreasonable taking into account that the B-factories 
have measured double charm and charmonium cross- 
sections. For example, for the very exclusive channel 
e + e~ —> J/ip + rjc, Belle finds a cross section of about 
26(6) femtobarn, while Babar reports some 18(5) femto- 
barn. Exclusive triple charmonium channels should be 
another two to three orders of magnitude smaller, but 
open flavor channels as the one we propose in Eq. (46) 
will be affected less by wave- function suppression, if the 
accelerator reaches sufficient energy. 

Adding the masses of all the particles, we find a 
threshold 

Mn acc + M p + 3il% ~ 11440(250) MeV (47) 

which lies at a slightly higher energy than Belle's data 
base at the T(5S), taken around 10860 MeV. 

It should be feasible for Super-B and Belle-II to take 
data at slightly higher energies around the T(6S) reso- 
nance and try to identify a triply charmed spectrum. 

As for the detection of an effect of the three-body 
force in the ground state spectrum, we do not share 
the optimism in Ref. [38], since the 25 MeV effect has 
to be found by comparing experimental data with de- 
tailed theoretical predictions of the masses, that, as we 
have shown, have systematic errors larger by an order 
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of magnitude. A cleaner way of extracting this informa- 
tion will have to be devised. 



11 Outlook 

We believe that we have presented a stride beyond the 
computation of Jia for triply heavy baryons. In this ini- 
tial evaluation we have employed the newly available 
NLO and NNLO potentials, obtaining results broadly 
consistent with other approaches, for the static poten- 
tial, and also with recoil corrections. 

Although at this point our prediction of the fl ccc 
mass at 4900(250) MeV is not particularly precise, it is 
information stemming from pNRQCD, the appropriate 
effective theory of QCD for ground state triply-heavy 
baryons, a qualitative improvement over the present, 
largely model-based situation. Computations by lattice 
groups are also underway. 

In future work we intend to reduce the uncertain- 
ties in this work by employing an additional, more so- 
phisticated Renormalon Subtracted scheme (RS) and 
by attempting a multi-wavefunction systematic diago- 
nalization. These two undertakings should address the 
biggest sources of uncertainty in the present work, the 
contribution to the binding of the infrared region, and 
the variational approximation. 

We have also been able to estimate the perturbative 
three-body force to be small. This is crucial information 
for the Faddeev formulation of the three-quark problem, 
since three-body forces are totally neglected there. At 
least we can now state that the error for the ground 
state, in the limit in which all three quarks are heavy, 
is modest and of order 25 MeV in the mass. 
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A Static potential with BLM scale fixing 

In this section we introduce the Brodsky-Lepage-McKenzie 
scale fixing at NLO (that we have used in practical computa- 
tions) as well as a sketch on how to proceed at NNLO (that 
we, however, have not further pursued). 



A.l NLO potential 

The idea of the BLM renormalization scale choice [58] is to ab- 
sorb the non-conformal terms proportional to the /3-function 
of QCD into the running coupling. The resulting scale /^blm 
is therefore fixed at NLO by demanding that the /3q term 
cancels. 

To one-loop order the idea can easily be carried out by 
shifting the scale of the strong coupling constant. This is done 
by substituting the solution of the renormalization group 
equation (with /3 ~ /3q taken at one-loop) 



l r l 2 1 



a " (A'blm) 



A>iog(^^ 



(48) 



«s(f*BL M ) + 



c<s (/'blm)' 



-in 



A> log(Mi 



r|r| 2 ) + . 



Demanding therefore that 
/3 (log( M | LM |r| 2 ) + 2 7B ) = 
we find that 



Mblm 



0.56- 



(49) 



For charmonium, one would estimate numerically that 



he 



0.75 GeV 



with the BLM scale correspondingly smaller, about 400 MeV 
as is usual. 

In terms of this scale, the potential to NLO takes the 
simpler form 



V LO + V NLO = [l + ai ) 



(50) 



A.2 NNLO potential 

The NNLO contribution to the static potential however con- 
tinues depending on /3q and fi\ 



V 



(0) 



■ ai- 



4?r 



(51) 



To eliminate these dependencies one needs to fix the BLM 
scale at two loops. However the expansion of the running 
coupling at the next order is significantly more difficult than 
Eq. (48). Then, to two loops, one has 



^ NLO (Q 2 ) 



(Q*)-b'a»^(Q*r log ( log 



A 2 
(52) 



where 



b' = 



Pi 



153 - lW f 



47T/3 2tt(33 - 2N f ) 

Since the double logarithm will yield a transcendental 
equation, one would have to fix the scale numerically. The 
computer algorithm proceeds as follows 



Triply heavy baryons in pNRQCD 



19 



1. Fix the BLM scale analytically at NLO as in Eq. (49) 
above. Then ^blm and LO {y? BLM ) are given. 

2. Obtain A from Eq. (19). 

3. Obtain Mb£m° 2 from E< 3- ( 52 ) h V substituting it for Q 2 
there. The equation relates iU^lm with [tg^M and the 
dependency has to be solved for at the same time that 
one tries to make vanish all /9q and fii -proportional terms 
in the potential. This is best performed by an iterative 
Newton's method. 

We have found this unpractical for the time being and have 
only pursued the BLM method at NLO. 



B Numerical methods 

B.l Fourier transform 

To numerically transform potentials between momentum and 
coordinate space we employ a standard Fast Fourier trans- 
form algorithm that implements the discrete formula 



After reducing the Hamiltonian to a numerical matrix, this 
is diagonalized. Since the radial problem is one-dimensional, 
accuracy in the diagonalization is not an issue. 

In the LO evaluation of the potential the coupling con- 



stant is fixed at a renormalization scale fj, 



or fj, 



employing the NLO running. Except for this small modifica- 
tion (needed since the same coupling constant cannot sensibly 
be used in both charmonium and bottomonium systems), the 
computation is perfectly consistent with perturbation theory, 
so that at NNLO, in the pure NNLO potential pieces the 
LO coupling constant is employed, whereas in the LO piece 
the NNLO coupling constant features, and so forth. As far as 
we can imagine no contamination is introduced from higher 
pieces in perturbation theory. We have also employed the per- 
turbative formulae for rj/jp-nj . 

We have of course checked the sensitivity of the numer- 
ical results to the number of points used in the grid solving 
Schroedinger's equation (300 turns out to yield very precise 
answers for low-lying states in the respective potentials), the 
maximum size of the grid (that extends to 4 fm and beyond) 
and other numeric artifacts. 



i)0" -i)M 



V N 



The wanted continuous transform is 



V(r) = 



d 3 q 
(2tt) 3 



? V(|q|) 



(53) 



(54) 



that after performing the angular integrals and grouping terms 
becomes 



V(r) = 



(2*-) = 



-Re 



i I qdqV(q)e tqr 



(55) 



We discretize a momentum interval (e, A) where A ~ 
50 GeV is well above the quark mass scale (the hard scale) 
and e is of order (20 fm) ~ 1 , well below any soft scale in the 
problems treated. The momentum variable is then stepped 
linearly according to q-j = A ~" j. The conjugate coordinate 
variable is automatically discretized as r n = 2nn/(A — e). 

With these choices the vectors that appear in Eq. (53) 

are 



V A An) 



{A- 



N 2 
-1 



-jexp 



2ttj 
If 



0-1) )V(qj) 



(56) 



2n 2 ri 



eX P( Jf l > Al 



B.2 Minimization 

We have written a computer programme that employs the 
well-known Minuit minimization package from CERN [59] 
to fix the values of a a , m c , and mj with the best possi- 
ble description of the observables that we have selected. The 
Schroedinger equation for the reduced particle is solved quasi- 
exactly (on a computer) with the perturbative potential to 
LO, NLO and NNLO, for both charm and bottom quarks. 
This is performed by discretizing the second derivative of the 
reduced radial function with the symmetric formula 



'(n) ~ 



u(r i+1 ) + u(fi-i) - 2u(r-j) 
h 2 



B.3 Iterative scale determination 

When the argument of the coupling constant depends on the 
coupling constant itself such as ot s = St s (ma s ) an iterative 
method is in order. 

We employ Newton's iterative numerical method. Denot- 
ing as as the successive approximations to the numeric 
value, and a s the function in Eq. (52) or similar, then the 
coupling constant has been found when 



F(a s ) 



M 



J J ill 



a s — 



(57) 



Newton's iteration, as long as F with sufficient sig- 
nificance, is given by 



F(a 



(n-l)x 



F'{a 



') 



(58) 



One can take as initial guess a° = & a (fj, 2 ) at any standard 
renormalization scale, and then iterate Eq. (58). 

As an alternative we also employ Jacobi's fixed point 



method, in which one starts with a guess ai°^ (presumably) 
larger than the true value, and then iterate the recursive re- 

until convergence. 



lation ai = a 



Mj,*, (0) 



B.4 Montecarlo computation of three-body matrix 
elements 

The variational matrix elements in the baryon computation 
are multidimensional integrals. Three particles, after center 
of mass separation, require six momentum integrations in the 
kinetic energy evaluation. Two-body potentials add one loop 
to the matrix element, up to nine dimensions. In the three- 
body force computation there are two exchanged momenta, 
and thus twelve-dimensional integrals. We make no attempt 
at separating rigid rotations and evaluate all these matrix 
elements numerically employing the Vegas algorithm [60,61]. 
In the matrix element 



(g) 
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we compute both numerator and denominator (the wavefunc- 
tion normalization) numerically. This makes the limits of inte- 
gration quite irrelevant for the computation since the function 
is normalized to one in the same region where the Hamilto- 
nian matrix element is computed, so that no probability den- 
sity is missed. In practice we never extend integration beyond 
the hard scale m c or rrif, (one expects the variational param- 
eters OL p and a\ to concentrate the momentum wavefunction 
around the soft scale ct s m c or a a mt). 

We employ a minimum of seven million evaluations of 
the Hamiltonian and reach a precision of about 10 MeV for 
standard computations, increasing this as needed. A 3 GHz 
processor can swipe an 8 x 8 set of variational parameters a p , 
a\ in about an hour. 

Our program performs wavefunction symmetrization (or 
mixed symmetrization for the ccb and bbc systems) by in- 
voking the wavefunction with exchanged spins and momen- 
tum arguments as needed. Although here we have not taken 
spin corrections into account, since they are unknown for the 
three-body problem, our program is performing (trivial) spin 
sums to allow for a simple upgrade once the spin kernels for 
triply heavy baryons become available. 
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